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Abstract 

English version: Because of their capability to preserve steady-states, 
well-balanced schemes for Shallow Water equations are becoming popular. 
Among them, the hydrostatic reconstruction proposed in [1] , coupled with 
a positive numerical flux, allows to verify important mathematical and 
physical properties like the positivity of the water height and, thus, to 
avoid unstabilities when dealing with dry zones. In this note, we prove 
that this method exhibits an abnormal behavior for some combinations of 
slope, mesh size and water height. 

Version Frangaise : Une limitation de la reconstruction hydro- 
statique pour la resolution du systeme de Saint- Venant. De 
par leur capacite a preserver les etats d'equilibre, les schemas equilibres 
connaissent actuellement un fort developpement dans la resolution des 
equations de Saint- Venant. En particulier, la reconstruction hydrosta- 
tique proposee dans [1], couplee a un flux numerique positif, permet de 
garantir certaines proprietes comme la positivite de la hauteur d'eau et, 
done, d'eviter certaines instabilites pour traiter les zones seches. Dans 
cette note, nous montrons que cette methode presente un defaut pour 
certaines combinaisons de pente, taille de maillage et hauteur d'eau. 



Version frangaise abregee 

Les equations de Saint- Venant (1) posent des difncultcs numeriqucs specifiqucs : 
preservation des equilibres stationnaires (flaques d'eau, lacs) et de la positivite 
de la hauteur d'eau. La reconstruction hydrostatique, introduite dans [1, 4], 
s'est imposee comme une methode particulierement adaptee. Elle fait partie de 
la classe des schemas dits equilibres (well-balanced). Partant d'une methode 
de volumes finis (2), dont le flux numerique F est adapte au systeme sans to- 
pographie, on reconstruit les variables u, h, h + z afin de preserver les equilibres. 
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Cette methode est appliquee le cas echeant a des variables deja reconstruites afin 
d'augmenter l'ordre de convergence. La mise en ceuvre complete de ce schema 
dans le cas d'vme reconstruction lineaire de type MUSCL est donnee par (3)-(7). 

Une reformulation de la reconstruction hydrostatique (8) met en evidence un 
comportement anormal pour certaines combinaisons de pente, maillage et hau- 
teur d'eau. Plus prcciscmcnt, il est montre que le schema surestime la hauteur 
d'eau dans les regions ou la relation (9) est vcrificc. 

Ce defaut est illustre par une serie de tests numcriqucs sur une solution 
analytiquc (voir 3). II s'avere particulicrcment spectaculaire a l'ordre 1, oil la 
methode calcule une meme pente apparente pour diffcrcntes valours effectives 
(Figure 2(Left)), mais il reste observable a l'ordre 2 (Figure 2(Right)). 

1 Introduction 

Derived from Navier-Stokes equations, Shallow Water equations describe the 
water flow properties as follows 

d t h + d x (hu) = 0, dt (hu) + d x (hu 2 + gh 2 /2) = -ghd x z, (1) 

where the unknowns are the water height h(t,x), the velocity u(t,x), and the 
topography z(x) is a given function. In what follows, we note the discharge 
q = hu, the vector of conservative variables U = (h hu)* and the flux F(U) = 
(hu gh 2 /2 + hu 2 y. The steady state of a lake at rest, or a puddle, (h + z = Cst 
and u = q = 0) is a particular solution to (1). Since [2], it is well known that 
the topography source term needs a special treatment in order to preserve at 
least this equilibrium. Such schemes are said to be well-balanced (since [9]). 

In the following, we present briefly the so-called hydrostatic reconstruction 
method, which permits, when coupled to a positive numerical flux, to obtain a 
family of well-balanced schemes that can preserve the water height nonncgativity 
and deal with dry zones. We show that this method, presented in [1, 4] and 
widely used, fails for some combinations of slope, mesh size and water height. 
We give the criteria that ensures the accuracy of the results. 

2 The numerical method 

The hydrostatic reconstruction follows the general principle of reconstruction 
methods. We start from a first order finite volume scheme for the homogeneous 
form of system (1): choosing a numerical flux F(Ul,Ur) (e.g. Rusanov, HLL, 
VFRoe-ncv, kinetic), a finite volume scheme takes the general form 

U* = U? - ^ [F(U i: U i+1 ) - P(Ui_ l5 Ui)] , (2) 

where At is the time step and Aa; the space step. Now the idea is to modify 
this scheme by applying the flux to reconstructed variables. Reconstruction 
can be used to get higher order schemes, in that case higher order in time 
is achieved through TVD-Rungc-Kutta methods. The aim of the hydrostatic 
reconstruction, which is described in the next section, is to be well-balanced, in 
the sense that it is designed to preserve at least steady states at rest (u = 0). 
When directly applied on the initial scheme, it leads to order one methods, while 
coupling it with high order accuracy reconstruction increases the order. 
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2.1 The hydrostatic reconstruction 

Wc describe now the implementation of this method for high order accuracy. 
The first step is to perform a high order reconstruction (MUSCL, ENO, UNO, 
WENO). To deal properly with the topography source term d x z, this recon- 
struction has to be performed on u, h and h + z, see [4]. This gives us the 
reconstructed values (£/_, Z-) and (£/+, z+), on which we apply the hydrostatic 
reconstruction [1, 4] on the water height, namely 

K+1/2L = max(/i i+1/2 _ + Zi+i/2- - max(z i+1/2 _, z i+1/2+ ), 0), 

Ui+i/2L = \hi+\/2Li ^i+l/2L u i+l/2-)> /o\ 

hi+i/2R = max(/i l+1/2+ + ^+i/2+ - max(z< + i/ 2 -,2i+i/2+), 0), 

, Ui+1/2R = (hi+l/2R,h i+ i/2RU i+ i/ 2 +)- 

With the hydrostatic reconstruction, the finite volume scheme (2) is modified 
as follows 

U; = U? At<S>(U n ) = U? ^ [F? +1/2L - FP_ 1/2R - Fc?] , (4) 



where 



r i+l/2L — r i+l/2 "+" J i+1/2L' r i-l/2R ~ r i-l/2 "+" °i-l/2R 1°; 



are left (respectively right) modifications of the initial numerical flux for the 
homogeneous problem. In this formula the flux is now applied to reconstructed 
variables: FV^_ l , 2 = F(C/ J '^ 1 ^ 2L , ^JYi/2i?,)' an< ^ we nave introduced 

/ \ / \ 

Sl+1/2L ~ \ f W +1/2 _ h1 +1/2L ) ) ' ^-V™ - ^ |(^_ 1/2+ - ^ 1/2R ) J 

(6) 

Finally, a centered source term has to be added to preserve consistency and 
well-balancing (see [1, 4]): 

( ° \ 

Fci = hi_ 1/2+ + /ij+i/2- , > ■ (7) 

V -3 2 (^+1/2- - ^-1/2+) J 

Formula (7) preserves the second order accuracy, but has to be modified for 
higher order approximations. 



2.2 The hydrostatic reconstruction rewritten 

We define Az i+1 / 2 = ^+1/2+ — Zj+i/2-- Once the space step and the high order 
reconstruction chosen, this is a fixed sequence. Now, the reconstructed variables 
(3) write 

f K+1/2L = max(/i i+ i/ 2 _ - max(0, Az i+1/2 ), 0), , g . 
\ h i+1/2 R = max(ft i+ i/ 2+ + min(0, Az l+1/2 ), 0). 

With (8), the defect of the hydrostatic reconstruction becomes apparent. To fix 
the ideas, suppose the local slope is positive, hence Az i+1 / 2 > 0, as in Figure 1. 
Then, for all Azj +1 / 2 such that Az i+1 / 2 > ^i+1/2- > 0; the reconstruction gives 
^i+i/2L — 0, while h i+ i/2u remains unchanged. In that case, the reconstruction 
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Figure 1: Default of the hydrostatic reconstruction. Left: threshold non acti- 
vated (limit case). Right: threshold activated 



prevents an unphysical negative value for h i+1 / 2 L, the counterpart of this being 
an underestimated difference h i+ i/ 2 L — ^i+i/afl- Therefore, there is a lack in 
the numerical flux computed from the modified Ricmann problem, which gives 
an underestimated velocity and consequently an overestimated height. This can 
be interpreted in terms of reconstructed slope as well, which is underestimated. 

In the general case, we can write the following local criterion of "non validity" 
for the hydrostatic reconstruction. 

Proposition 2.1 For a fixed discretization, if for some io < i < ii one has 

Azt+i/2 > K+i/2- > 0, or - Az l+1/2 > K +l/2+ > 0, (9) 

then the hydrostatic reconstruction will overestimate (resp. underestimate) the 
water height (resp. velocity). 

Notice that from a theoretical viewpoint, this is not limiting. Indeed, since 
this class of schemes is consistent with the system of partial differential equations 
(see [4]), the problem disappears when refining the discretization. However, it 
has to be taken into account for practical computations, with a fixed discretiza- 
tion. It is particularly apparent for order 1 schemes, but also remains present 
for order 2. 



3 Numerical illustration 

For numerical illustration purpose, we introduce an explicit solution which con- 
sists in a supercritical steady flow on an inclined plane with constant slope 
d x z = a (it is referenced in the SWASHES library [8]). Steady states are solu- 
tions to 

q(x, t)=q = Cst, d x (hu 2 + 9~^j = ~ghd x z, 

and the height profile h(x) must be a solution to Bernoulli's law rewritten as a 
third order equation in h: 

a3+ " ! (—4-' , °) + I^' (io) 

where ho and qo = hotio are the height and discharge at x = 0, which completely 
determine the profile since the flow is supercritical. A careful study of the roots 
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of the polynomial shows that the supercritical height profile is decreasing in x: 
h(x) < ho for all x > 0. 

The numerical strategy we choose consists in the HLL flux and a modified 
MUSCL reconstruction. In [7], this combination of flux and linear reconstruction 
has shown to be the best compromise between accuracy, stability and CPU 
time cost. We refer to [4, 7] for a presentation of the HLL flux. The MUSCL 
reconstruction of a scalar function s £ R is 

Ax Ax 

s i-l/2+ — s i 2~ s i+l/2- — s i H DSi, (11) 

where the "minmod" operator D is given by 

/ - ■ ■ - A fmin(x,j/) ifx,y>0, 

Dsi = minmod ( 1 — 1 — , J , minmod(x, y) = < max(x, y) if x, y < 0, 

[o else. 

(12) 

In order to keep the discharge conservation, the reconstruction of the velocity 
has to be modified as 

h i+ i/ 2 - Ax hi-U2+ Ax 
Ui-i/2+ = ut — — -Dm u l+1/2 - = Ui H — —Dm. 

Notice that if we take Dsi = in (11), then z i+1 / 2 - = ^i+i/2- = so that the 
centered term (7) disappears, and we recover the first order scheme in space. 
Second order in time is achieved through a classical Hcun predictor-corrector 
method. 

We turn now to the specific data for simulations. The domain length is 
L = 10 m, and we choose initial data 

ho = h(x = 0,t) = 0.02 m, q = q{x = 0,t) = 0.01 m 2 /s, 

which is indeed supercritical. All simulations are performed with a space step 
Air = 0.1 m, the time step is fixed in order to satisfy the CFL condition. 

The analytical solution is computed with 7 negative slopes a = 5%, 13%, 
14%, 15%, 16%, 17% and 18%, both at first and second order. Numerical results 
are compared to the analytical results in Figure 2, where a part of the domain 
is displayed (x between 1 m and 3 m). 

With this set of data, a domain of admissible slopes can be estimated for the 
order 1 hydrostatic reconstruction. Indeed since Az J+1 / 2 = |a|Ax, inserting a 
characteristic height of the flow h* in (9) gives a bound for the slope. One can 
use h* = ho, the incoming height. Since the height profile is decreasing, the 
whole profile will be wrong if ho < |a|Ax, that is here \a\ > 20%. But actually, 
since the height decreases quite rapidly, a more accurate estimate is obtained 
for x = 1.5 m with h* = 0.005 m (see Fig. 2), which leads to slopes above 5%. 

We observe on Figure 2 (Left) that with the first order scheme, the effect of 
the hydrostatic reconstruction is so important that all curves for slopes between 
13% and 18% are superposed. In that case, the apparent result is the simulation 
of a single slope, namely 13%. For a 5% slope, we still observe a slightly overes- 
timated water height, as anticipated since 5% is a limit case as observed above. 
With the second order, the water heights are still overestimated, but in a very 
slighter way, and the different curves are no longer identical (Figure 2 (Right)). 
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Figure 2: Default of hydrostatic reconstruction: water height at first and second 
order of accuracy for different slopes. Dotted curves are simulations, plain curves 
are exact solutions. Left: first and second order curves. Right: zoom on second 
order curves 

4 Conclusions 

The hydrostatic reconstruction may fail for certain combinations of water height, 
slope and mesh size, namely in regions where (9) holds. The defaults are par- 
ticularly apparent for order 1 schemes, leading to wrongly estimated slopes, but 
still remain at order 2, with some overestimated water heights. We emphasize 
that the problem disappears when refining the mesh, but has to be taken into 
account for a given discretization. The generalization of the hydrostatic recon- 
struction proposed in [6] does not exhibit the limitation discussed here, even 
for first order schemes, but positivity is not ensured. Other schemes involving 
threshold values (e.g. [5, 10]) very likely encounter the same kind of problem. 
Alternatively, the scheme recently introduced in [3] preserves the water height 
positivity and does not suffer from this problem. Notice finally that criterion 
(9) may be of some utility for adaptive mesh schemes, such as the ones used in 
Gerris [11]. 
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